#ifndef __SPLINE__
#define __SPLINE__

#include <iostream>
#include <cmath>
#include <algorithm>
#include <vector>
#include <eigen3/Eigen/Dense>
#include <eigen3/Eigen/Sparse>
#include "Function.h"
#include "Polynomial.h"
#include "Interpolation.h"

using namespace std;
using namespace Eigen;
class ppForm_interpolation : public Interpolation
{
private:
  Function &f;
  vector<double> Point;
  vector<Polynomial> Poly;
  vector<vector<double>> Coef;
  int n, Boundary_conditions, order;
  VectorXd m;
public:
  ppForm_interpolation(Function &_f, vector<double> _Point, int _Boundary_conditions,int _order): f(_f), Point(_Point), Boundary_conditions(_Boundary_conditions), order(_order)
  {
    if ((Boundary_conditions != 1) && (Boundary_conditions != 2) && (Boundary_conditions != 3) )
      {
	cout<< "error: There is no such boundary condition." <<endl;
	  exit(0);
      }
    if ((order != 1) && (order != 3))
      {
	cout<< "error: There is no such order." <<endl;
	exit(0);
      }
    n = Point.size();
  }

  double operator()(double _x)
  {
    if ((_x < Point[0]) || (_x > Point[n-1]))
      {
	cout<< "error: out of range." <<endl;
	exit(0);
      }
    int i = 0;
    while (_x > Point[i])
      {
	i++;
      }
    if (_x == Point[i])
      {
	return f(Point[i]);
      }
    else
      {
	return Poly[i-1](_x);
      }
  }

  void solve()
  {
    if (order == 1)
      {
	vector<double> C;
	for (int i = 1; i < n; i++)
	  {
	    C.clear();
	    C.push_back(f(Point[i]) - Point[i] * (f(Point[i]) - f(Point[i-1]))/(Point[i] - Point[i-1]));
	    C.push_back((f(Point[i]) - f(Point[i-1]))/(Point[i] - Point[i-1]));
	    Coef.push_back(C);
	    Polynomial p(C,0);
	    Poly.push_back(p);
	  }
	
      }
    else
      {
	SparseMatrix<double> A(n,n);
	MatrixXd b(n,1);
	vector<Triplet<double>> triplets;
	double y,z;

	b = MatrixXd::Zero(n,1);
	for (int i = 1;i < n-1; i++)
	  {
	    y = (Point[i+1] - Point[i])/(Point[i+1] - Point[i-1]);
	    z = (Point[i] - Point[i-1])/(Point[i+1] - Point[i-1]);
	    triplets.push_back(Triplet<double>(i, i, 2.0));
	    triplets.push_back(Triplet<double>(i, i-1, y));
	    triplets.push_back(Triplet<double>(i, i+1, z));
	    b(i,0)=3.0*z*(f(Point[i+1])-f(Point[i]))/(Point[i+1]-Point[i]) + 3*y*(f(Point[i])-f(Point[i-1]))/(Point[i]-Point[i-1]);
	  }

	if (Boundary_conditions == 1 )
	  {
	    triplets.push_back(Triplet<double>(0, 0, 1.0));
	    triplets.push_back(Triplet<double>(n-1, n-1, 1.0));
	    b(0,0) = f.diff(Point[0]);
	    b(n-1,0) = f.diff(Point[n-1]);
	  }
	else if (Boundary_conditions == 2)
	  {
	    triplets.push_back(Triplet<double>(0, 0, 4.0));
	    triplets.push_back(Triplet<double>(0, 1, 2.0));
	    triplets.push_back(Triplet<double>(n-1, n-1, 4.0));
	    triplets.push_back(Triplet<double>(n-1, n-2, 2.0));
	    b(0,0) = 6.0*(f(Point[1])-f(Point[0]))/(Point[1]-Point[0]) - f.diff2(Point[0])*(Point[1]-Point[0]);
	    b(n-1,0) = 6.0*(f(Point[n-1])-f(Point[n-2]))/(Point[n-1]-Point[n-2]) + f.diff2(Point[n-1])*(Point[n-1]-Point[n-2]);
	  }
	else
	  {
	    triplets.push_back(Triplet<double>(0, 0, 4.0));
	    triplets.push_back(Triplet<double>(0, 1, 2.0));
	    triplets.push_back(Triplet<double>(n-1, n-1, 4.0));
	    triplets.push_back(Triplet<double>(n-1, n-2, 2.0));
	    b(0,0) = 6.0*(f(Point[1])-f(Point[0]))/(Point[1]-Point[0]);
	    b(n-1,0) = 6.0*(f(Point[n-1])-f(Point[n-2]))/(Point[n-1]-Point[n-2]);
	  }

	A.setFromTriplets(triplets.begin(), triplets.end());
	A.makeCompressed();
	SparseLU<SparseMatrix<double> > Solver_sparse;
	Solver_sparse.compute(A);
	m = Solver_sparse.solve(b);

	for (int i = 1; i < n; i++)
	  {
	    double K;
	    vector<double> C;
	    K = (f(Point[i])-f(Point[i-1]))/(Point[i]-Point[i-1]);
	    C.push_back(f(Point[i-1]));
	    C.push_back(m(i-1));
	    C.push_back((3.0*K-2.0*m(i-1) - m(i))/(Point[i]-Point[i-1]));
	    C.push_back((m(i-1)+m(i)-2.0*K)/pow(Point[i]-Point[i-1],2));
	    Coef.push_back(C);
	    Polynomial p(C, Point[i-1]);
	    Poly.push_back(p);
	  }
      }
  }
};


class Bspline_interpolation : public Interpolation
{
private:
  Function &f;
  vector<double> Point,add1,add2;
  vector<B_spline> B;
  //vector<Polynomial> Poly;
  //vector<vector<double>> Coef;
  int n, Boundary_conditions, order;
  VectorXd m;
public:
  Bspline_interpolation(Function &_f, vector<double> _Point, int _Boundary_conditions,int _order): f(_f), Point(_Point), Boundary_conditions(_Boundary_conditions), order(_order)
  {
    if ((Boundary_conditions != 1) && (Boundary_conditions != 2) && (Boundary_conditions != 3) )
      {
	cout<< "error: There is no such boundary condition." <<endl;
	  exit(0);
      }
    if ((order != 1) && (order != 3))
      {
	cout<< "error: There is no such order." <<endl;
	exit(0);
      }
    n = Point.size();   
  }

  void solve()
  {
    if (order == 1)
      {
	add1 = {2 * Point[0] - Point[1]};
	add2 = {Point[n-1] + Point[1] - Point[0]};
      }
    else
      {
	add1 = {3 * Point[1] - 2 * Point[0], Point[0] -2 * Point[1], 2 * Point[0] - Point[1]};
	add2 = {Point[n-1] + Point[1] - Point[0], Point[n-1] + 2 * (Point[1] - Point[0]), Point[n-1] + 3 * (Point[1] - Point[0])};
      }
    Point.insert(Point.begin(), add1.begin(), add1.end());
    Point.insert(Point.end(), add2.begin(), add2.end());
    for (int i = 2 - order; i < n + 1; i++)
      {
	vector<double> h(Point.begin() + i - 2 + order, Point.begin() + i + 2 * order);
	B_spline b(order, h);
	B.push_back(b);
      }

    if (order == 1)
      {
	SparseMatrix<double> A(n, n);
	vector<Triplet<double> > triplets;
	MatrixXd b(n, 1);
	b = MatrixXd::Zero(n, 1);
	for (int i = 0; i < n; i++)
	  {
	    triplets.push_back(Triplet<double> (i, i, 1.0));
	    b(i, 1) = f(Point[i + 1]);
	  }
	A.setFromTriplets(triplets.begin(), triplets.end());
	A.makeCompressed();
	SparseLU<SparseMatrix<double> > Solver_sparse;
	Solver_sparse.compute(A);
	m = Solver_sparse.solve(b); 
      }
    else
      {
	SparseMatrix<double> A(n+2, n+2);
	vector<Triplet<double> > triplets;
	MatrixXd b(n+2, 1);
	b = MatrixXd::Zero(n+2, 1);
	for (int i = 1; i < n; i++)
	  {
	    triplets.push_back(Triplet<double>(i, i-1, B[i-1](Point[i+2])));
	    triplets.push_back(Triplet<double>(i, i-1, B[i](Point[i+2])));
	    triplets.push_back(Triplet<double>(i, i-1, B[i+1](Point[i+2])));
	  }
	if (Boundary_conditions == 1)
	  {
	    b(0, 0) = f.diff(Point[3]);
	    b(n+1,0) = f.diff(Point[n+2]);
	    B_spline b0(2, vector<double> (Point.begin()+1, Point.begin()+5));
	    B_spline b1(2, vector<double> (Point.begin()+2, Point.begin()+6));
	    B_spline bn(2, vector<double> (Point.end()-5, Point.end()-1));
	    B_spline bn_1(2, vector<double> (Point.end()-6, Point.end()-2));
	    triplets.push_back(Triplet<double>(0, 0,(-3.0*b0(Point[3]))/(Point[4]-Point[1])));
	    triplets.push_back(Triplet<double>(0, 1,(3.0*b0(Point[3]))/(Point[4]-Point[1]) - (3.0*b1(Point[3]))/(Point[5]-Point[2])));
	    triplets.push_back(Triplet<double>(0, 1,(3.0*b1(Point[3]))/(Point[5]-Point[2])));
	    triplets.push_back(Triplet<double>(n+1, n-1,(-3.0*bn_1(Point[n+2]))/(Point[n+3]-Point[n])));
	    triplets.push_back(Triplet<double>(n+1, n,(3.0*bn_1(Point[n+2]))/(Point[n+3]-Point[n]) - (3.0*bn(Point[n+2]))/(Point[n+4]-Point[n+1])));
	    triplets.push_back(Triplet<double>(n+1, n+1,(3.0*bn(Point[n+2]))/(Point[n+4]-Point[n+1])));
	  }
	else if (Boundary_conditions == 2)
	  {
	    b(0, 0) = f.diff2(Point[3]);
	    b(n+1,0) = f.diff2(Point[n+2]);
	    triplets.push_back(Triplet<double>(0, 0, 6.0/((Point[4]-Point[1])*(Point[4]-Point[2]))));
	    triplets.push_back(Triplet<double>(0, 1, -6.0/((Point[4]-Point[1])*(Point[4]-Point[2])) - 6.0/((Point[5]-Point[2])*(Point[4]-Point[2]))));
	    triplets.push_back(Triplet<double>(0, 2, 6.0/((Point[5]-Point[2])*(Point[4]-Point[2]))));
	    triplets.push_back(Triplet<double>(n+1, n-1, 6.0/((Point[n+3]-Point[n])*(Point[n+3]-Point[n+1]))));
	    triplets.push_back(Triplet<double>(n+1, n, -6.0/((Point[n+3]-Point[n])*(Point[n+3]-Point[n+1])) - 6.0/((Point[n+4]-Point[n+1])*(Point[n+3]-Point[n+1]))));
	    triplets.push_back(Triplet<double>(n+1, n+1, 6.0/((Point[n+4]-Point[n+1])*(Point[n+3]-Point[n+1]))));
	  }
	else
	  {
	    b(0, 0) = 0.0;
	    b(n+1,0) = 0.0;
	    triplets.push_back(Triplet<double>(0, 0, 6.0/((Point[4]-Point[1])*(Point[4]-Point[2]))));
	    triplets.push_back(Triplet<double>(0, 1, -6.0/((Point[4]-Point[1])*(Point[4]-Point[2])) - 6.0/((Point[5]-Point[2])*(Point[4]-Point[2]))));
	    triplets.push_back(Triplet<double>(0, 2, 6.0/((Point[5]-Point[2])*(Point[4]-Point[2]))));
	    triplets.push_back(Triplet<double>(n+1, n-1, 6.0/((Point[n+3]-Point[n])*(Point[n+3]-Point[n+1]))));
	    triplets.push_back(Triplet<double>(n+1, n, -6.0/((Point[n+3]-Point[n])*(Point[n+3]-Point[n+1])) - 6.0/((Point[n+4]-Point[n+1])*(Point[n+3]-Point[n+1]))));
	    triplets.push_back(Triplet<double>(n+1, n+1, 6.0/((Point[n+4]-Point[n+1])*(Point[n+3]-Point[n+1]))));
	  }
	A.setFromTriplets(triplets.begin(), triplets.end());
	A.makeCompressed();
	SparseLU<SparseMatrix<double> > Solver_sparse;
	Solver_sparse.compute(A);
	m = Solver_sparse.solve(b); 
      }
  }

  double operator()(double _x)
  {
    if (order == 1)
      {
	if ((_x < Point[1]) || (_x > Point[n]))
	  {
	    cout<< "error: out of range" <<endl;
	    exit(0);
	  }
	if (_x == Point[1])
	  {
	    return f(Point[1]);
	  }
	int j = 1;
	while (_x > Point[j+1])
	  {
	    j++;
	  }
	return m(j-1)*B[j-1](_x) + m(j)*B[j](_x);
      }
    else
      {
	if ((_x < Point[3]) || (_x > Point[n+2]))
	  {
	    cout<< "error: out of range" <<endl;
	    exit(0);
	  }
	if (_x == Point[3])
	  {
	    return f(Point[3]);
	  }
	int j = 1;
	while (_x > Point[j+1])
	  {
	    j++;
	  }
	return m(j-1)*B[j-1](_x) + m(j)*B[j](_x) + m(j+1)*B[j+1](_x) + m(j+2) * B[j+2](_x);
      }
    return 0;
  }
  
};
  
#endif
    
  

  
